!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

!==================================================== John C. Warner ===
!                                                                      !
!  This routine computes floc transformation.                          !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Verney, R., Lafite, R., Claude Brun-Cottan, J., & Le Hir, P. (2011).!
!    Behaviour of a floc population during a tidal cycle: laboratory   !
!    experiments and numerical modelling. Continental Shelf Research,  !
!    31(10), S64-S83.                                                  !
!=======================================================================

Module Mod_FlocMod
#if defined (SEDIMENT) && (CSTMS_SED)

# if defined (WAVE_CURRENT_INTERACTION)
USE MOD_WAVE_CURRENT_INTERACTION
# endif

Use All_vars
Use Mod_Par
Use Mod_Prec 
Use Mod_Types


Use Mod_CSTMS_vars

! 
!  switched to activate the flocmod model
!
  implicit none
  
  public

!  l_ADS           : boolean set to .true. if differential settling aggregation 
!
!  l_ASH           : boolean set to .true. if shear aggregation
!
!  l_COLLFRAG      : boolean set to .true. if collision-induced fragmentation enable
!  f_dp0           : Primary particle size (m), typically 4e-6 m
!
!  f_nf            : Floc fractal dimension, typically ranging from 1.6 to 2.6
!
!  f_dmax          : (unused) Maximum diameter (m)
!
!  f_nb_frag       : number of fragments by shear erosion. If binary/ternary : 2.
!
!  f_alpha         :  flocculation efficiency, ranging from 0 .to 1.
!
!  f_beta          : shear fragmentation rate
!
!  f_ater          : for ternary breakup, use 0.5, for binary : 0.
!                   (a boolean could be better)
!
!  f_ero_frac      : fraction of the shear fragmentation term 
!                    transfered to shear erosion. 
!                    Ranging from 0. (no erosion) to 1. (all erosion)
!
!  f_ero_nbfrag    : Number of fragments induced by shear erosion.   
!
!  f_ero_iv        : fragment size class (could be changed to a particle 
!                    size or a particle distribution (INTEGER)
!
!  f_collfragparam : Fragmentation rate for collision-induced breakup
!
! f_clim           : min concentration below which flocculation
!                    processes are not calculated
!
! l_testcase - if .TRUE. sets G(t) to values from Verney et al., 2011 lab experiment
!
! GLS_P    : Stability exponent (non-dimensional).
!
! GLS_M    : Turbulent kinetic energy exponent (non-dimensional).
!
! GLS_N    : Turbulent length scale exponent (non-dimensional).
!
! GLS_CMU0 : Stability coefficient.


  logical :: l_ASH
  logical :: l_ADS
  logical :: l_COLLFRAG
  logical :: l_testcase
  INTEGER  :: f_ero_iv
  real(sp) :: f_dp0,f_alpha,f_beta,f_nb_frag
  real(sp) :: f_dmax,f_ater,f_clim
  real(sp) :: f_ero_frac,f_ero_nbfrag
  real(sp) :: f_nf
  real(sp) :: f_frag
  real(sp) :: f_fter
  real(sp) :: f_collfragparam
  real(sp) :: gls_p
  real(sp) :: gls_m
  real(sp) :: gls_n
  real(sp) :: gls_cmu0
  real(sp), parameter :: rhoref = 1030.0_sp

  
  ! Dynamics selecting switches
  logical :: FLOC_TURB_DISS
  logical :: FLOC_BBL_DISS


  contains


  SUBROUTINE allocate_sedflocs

!  Imported variable declarations.
!
!
!-----------------------------------------------------------------------
!  Allocate structure variables.
!-----------------------------------------------------------------------
!
!
  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: allocate_sedflocs"

      allocate ( SEDFLOCS % f_diam(NCS) )
      allocate ( SEDFLOCS % f_vol(NCS) )
      allocate ( SEDFLOCS % f_rho(NCS) )
      allocate ( SEDFLOCS % f_cv(NCS) )
      allocate ( SEDFLOCS % f_l3(NCS) )
      allocate ( SEDFLOCS % f_mass(0:NCS+1) )
      allocate ( SEDFLOCS % f_coll_prob_sh(NCS,NCS) )
      allocate ( SEDFLOCS % f_coll_prob_ds(NCS,NCS) )
      allocate ( SEDFLOCS % f_l1_sh(NCS,NCS) )
      allocate ( SEDFLOCS % f_l1_ds(NCS,NCS) )
      allocate ( SEDFLOCS % f_g3(NCS,NCS) )
      allocate ( SEDFLOCS % f_l4(NCS,NCS) )
      allocate ( SEDFLOCS % f_g1_sh(NCS,NCS,NCS) )
      allocate ( SEDFLOCS % f_g1_ds(NCS,NCS,NCS) )
      allocate ( SEDFLOCS % f_g4(NCS,NCS,NCS) )
  if(dbg_set(dbg_sbr)) write(ipt,*) "End: allocate_sedflocs"

  RETURN

  END SUBROUTINE allocate_sedflocs


  SUBROUTINE initialize_sedflocs_param (f_mass,f_diam,f_g1_sh,          &
     &                              f_g1_ds,f_g3,f_l1_sh,f_l1_ds,       &
     &                              f_coll_prob_sh,f_coll_prob_ds,      &
     &                              f_l3)


  Use Scalar
  Use Control, only : ireport,iint,msr,par,pi
  Use Lims,    only : m,mt,nt,kbm1,numqbc,kb,nprocs,myid

  implicit none 
!
!  Imported variable declarations.
!
      real(sp), intent(inout) :: f_mass(0:NCS+1)
      real(sp), intent(inout) :: f_diam(NCS)
      real(sp), intent(inout) :: f_g1_sh(NCS,NCS,NCS)
      real(sp), intent(inout) :: f_g1_ds(NCS,NCS,NCS)
      real(sp), intent(inout) :: f_g3(NCS,NCS)
      real(sp), intent(inout) :: f_l1_sh(NCS,NCS)
      real(sp), intent(inout) :: f_l1_ds(NCS,NCS)
      real(sp), intent(inout) :: f_coll_prob_sh(NCS,NCS)
      real(sp), intent(inout) :: f_coll_prob_ds(NCS,NCS)
      real(sp), intent(inout) :: f_l3(NCS)
!
!  Local variable declarations.
!
   logical  :: f_test
   real(sp) :: f_weight,mult,dfragmax
   integer  :: iv1,iv2,iv3,iv,itrc
   real(sp) :: f_vol(NCS),f_rho(NCS)
   real(sp), parameter :: mu = 0.001_sp
   real(sp) :: eps 
   REAL(DP), PARAMETER :: g = 9.81_SP
   eps = epsilon(1.0)
!
! ALA the rest of the initialization
!
!      l_ADS=.false.
!      l_ASH=.true.
!      l_COLLFRAG=.false.
!      f_dp0=0.000004_sp
!      f_nf=1.9_sp
!      f_dmax=0.001500_sp
!      f_nb_frag=2.0_sp
!      f_alpha=0.35_sp
!      f_beta=0.15_sp
!      f_ater=0.0_sp
!      f_ero_frac=0.0_sp
!      f_ero_nbfrag=2.0_sp
!      f_ero_iv=1
!      f_collfragparam=0.01_sp

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: initialize_sedflocs_param"


!!--------------------------------------------------
!! floc characteristics
   DO itrc=1,NCS
      f_diam(itrc)=sed(itrc)%Sd50
      f_vol(itrc)=pi/6.0_sp*(f_diam(itrc))**3.0_sp
      f_rho(itrc)=rhoref+(2650.0_sp-rhoref)*                          &
    &     (f_dp0/f_diam(itrc))**(3.0_sp-f_nf)
      f_mass(itrc)=f_vol(itrc)*(f_rho(itrc)-rhoref)
   ENDDO
   f_mass(NCS+1)=f_mass(NCS)*2.0_sp+1.0_sp  
   IF (f_diam(1).eq.f_dp0)  THEN
       f_mass(1)=f_vol(1)*sed(1)%Srho
   ENDIF

   IF(MSR)THEN
!TODO - This wont parallelize 
      WRITE(*,*) ' '
      WRITE(*,*) '*** FLOCMOD INIT *** '
      write(*,*) 'NCS, NNS:',NCS,NNS   
      WRITE(*,*) 'class diameter (um)  volume (m3)  density (kg/m3)  mass (kg)'
      DO itrc=1,NCS
         WRITE(*,*) itrc,f_diam(itrc)*1e6,f_vol(itrc),f_rho(itrc),f_mass(itrc)
      ENDDO
      write(*,*) 'f_mass(0) and f_mass(NCS+1): ',f_mass(0),f_mass(NCS+1)
      WRITE(*,*) ' '
      WRITE(*,*) ' *** PARAMETERS ***'
      WRITE(*,*) 'Primary particle size (f_dp0)                                : ',f_dp0
      WRITE(*,*) 'Fractal dimension (f_nf)                                     : ',f_nf
      WRITE(*,*) 'Flocculation efficiency (f_alpha)                            : ',f_alpha
      WRITE(*,*) 'Floc break up parameter (f_beta)                             : ',f_beta
      WRITE(*,*) 'Nb of fragments (f_nb_frag)                                  : ',f_nb_frag
      WRITE(*,*) 'Ternary fragmentation (f_ater)                               : ',f_ater
      WRITE(*,*) 'Floc erosion (% of mass) (f_ero_frac)                        : ',f_ero_frac
      WRITE(*,*) 'Nb of fragments by erosion (f_ero_nbfrag)                    : ',f_ero_nbfrag
      WRITE(*,*) 'fragment class (f_ero_iv)                                    : ',f_ero_iv
!      WRITE(*,*) 'negative mass tolerated before redistribution (f_mneg_param) : ',f_mneg_param
      WRITE(*,*) 'Boolean for differential settling aggregation (L_ADS)        : ',l_ADS
      WRITE(*,*) 'Boolean for shear aggregation (L_ASH)                        : ',l_ASH
      WRITE(*,*) 'Boolean for collision fragmenation (L_COLLFRAG)              : ',l_COLLFRAG
      WRITE(*,*) 'Collision fragmentation parameter (f_collfragparam)          : ',f_collfragparam
      WRITE(*,*) ' '
      WRITE(*,*) 'Value of eps                                                 : ',eps
      WRITE(*,*) ' '
   END IF
! kernels computation former SUBROUTINE flocmod_kernels


!!--------------------------------------------------------------------------
!! * Executable part
      f_test=.true.
      dfragmax=0.00003_sp
! compute collision probability former SUBROUTINE flocmod_agregation_statistics

      DO iv1=1,NCS
        DO iv2=1,NCS

          f_coll_prob_sh(iv1,iv2)=1.0_sp/6.0_sp*(f_diam(iv1)+           &
     &            f_diam(iv2))**3.0_sp

          f_coll_prob_ds(iv1,iv2)=0.25_sp*pi*(f_diam(iv1)+              &
     &            f_diam(iv2))**2.0_sp*g/mu*abs((f_rho(iv1)-         &
     &            rhoref)*f_diam(iv1)**2.0_sp-(f_rho(iv2)-rhoref)*      &
     &            f_diam(iv2)**2.0_sp)

        ENDDO
      ENDDO

  !********************************************************************************
  ! agregation : GAIN : f_g1_sh and f_g1_ds
  !********************************************************************************
      DO iv1=1,NCS
       DO iv2=1,NCS
        DO iv3=iv2,NCS
           IF((f_mass(iv2)+f_mass(iv3)) .gt. f_mass(iv1-1)    .and.     &
     &            ((f_mass(iv2)+f_mass(iv3)) .le. f_mass(iv1))) THEN

              f_weight=(f_mass(iv2)+f_mass(iv3)-f_mass(iv1-1))/         &
     &            (f_mass(iv1)-f_mass(iv1-1))

           ELSEIF ((f_mass(iv2)+f_mass(iv3)) .gt. f_mass(iv1) .and.     &
     &            ((f_mass(iv2)+f_mass(iv3)) .lt. f_mass(iv1+1))) THEN

              IF (iv1 .eq. NCS) THEN
                 f_weight=1.0_sp
              ELSE
                 f_weight=1.0_sp-(f_mass(iv2)+f_mass(iv3)-              &
     &                 f_mass(iv1))/(f_mass(iv1+1)-f_mass(iv1))
              ENDIF

           ELSE
              f_weight=0.0_sp
           ENDIF

           f_g1_sh(iv2,iv3,iv1)=f_weight*f_alpha*                       &
     &            f_coll_prob_sh(iv2,iv3)*(f_mass(iv2)+                 &
     &            f_mass(iv3))/f_mass(iv1)
           f_g1_ds(iv2,iv3,iv1)=f_weight*f_alpha*                       &
     &            f_coll_prob_ds(iv2,iv3)*(f_mass(iv2)+                 &
     &            f_mass(iv3))/f_mass(iv1)

        ENDDO
       ENDDO
      ENDDO

  !********************************************************************************
  ! Shear fragmentation : GAIN : f_g3
  !********************************************************************************

      DO iv1=1,NCS
       DO iv2=iv1,NCS

        IF (f_diam(iv2).gt.dfragmax) THEN
           ! binary fragmentation

           IF (f_mass(iv2)/f_nb_frag .gt. f_mass(iv1-1) &
                .and. f_mass(iv2)/f_nb_frag .le. f_mass(iv1)) THEN

              IF (iv1 .eq. 1) THEN 
                 f_weight=1.0_sp
              ELSE
                 f_weight=(f_mass(iv2)/f_nb_frag-f_mass(iv1-1))/        &
     &                (f_mass(iv1)-f_mass(iv1-1))
              ENDIF

           ELSEIF (f_mass(iv2)/f_nb_frag .gt. f_mass(iv1)               &
     &            .and. f_mass(iv2)/f_nb_frag .lt. f_mass(iv1+1)) THEN

              f_weight=1.0_sp-(f_mass(iv2)/f_nb_frag-f_mass(iv1))/      &
     &              (f_mass(iv1+1)-f_mass(iv1))

           ELSE

              f_weight=0.0_sp

           ENDIF

        ELSE
           f_weight=0.0_sp
        ENDIF

        f_g3(iv2,iv1)=f_g3(iv2,iv1)+(1.0_sp-f_ero_frac)*(1.0_sp-        &
     &            f_ater)*f_weight*f_beta*f_diam(iv2)*((f_diam(iv2)-    &
     &            f_dp0)/f_dp0)**(3.0_sp-f_nf)*f_mass(iv2)/             &
     &            f_mass(iv1)

        ! ternary fragmentation
        IF (f_diam(iv2).gt.dfragmax) THEN
           IF (f_mass(iv2)/(2.0_sp*f_nb_frag) .gt. f_mass(iv1-1) .and.  &
     &          f_mass(iv2)/(2.0_sp*f_nb_frag) .le. f_mass(iv1)) THEN

              IF (iv1 .eq. 1) THEN 
                 f_weight=1.0_sp
              ELSE
                 f_weight=(f_mass(iv2)/(2.0_sp*f_nb_frag)-              &
     &                f_mass(iv1-1))/(f_mass(iv1)-f_mass(iv1-1))
              ENDIF

           ELSEIF (f_mass(iv2)/(2.0_sp*f_nb_frag) .gt. f_mass(iv1)      &
     &           .and. f_mass(iv2)/(2.0_sp*f_nb_frag) .lt.              &
     &            f_mass(iv1+1)) THEN

              f_weight=1.0_sp-(f_mass(iv2)/(2.0_sp*f_nb_frag)-          &
     &                 f_mass(iv1))/(f_mass(iv1+1)-f_mass(iv1))

           ELSE
              f_weight=0.0_sp

           ENDIF
           ! update for ternary fragments
           f_g3(iv2,iv1)=f_g3(iv2,iv1)+(1.0_sp-f_ero_frac)*(f_ater)*    &
     &            f_weight*f_beta*f_diam(iv2)*((f_diam(iv2)-f_dp0)/     &
     &            f_dp0)**(3.0_sp-f_nf)*f_mass(iv2)/f_mass(iv1)   

           ! Floc erosion

           IF ((f_mass(iv2)-f_mass(f_ero_iv)*f_ero_nbfrag) .gt.         &
     &            f_mass(f_ero_iv)) THEN

              IF (((f_mass(iv2)-f_mass(f_ero_iv)*f_ero_nbfrag) .gt.     &
     &            f_mass(iv1-1)) .and. (f_mass(iv2)-f_mass(f_ero_iv)*   &
     &            f_ero_nbfrag) .le. f_mass(iv1)) THEN

                 IF (iv1 .eq. 1) THEN
                    f_weight=1.0_sp
                 ELSE
                    f_weight=(f_mass(iv2)-f_mass(f_ero_iv)*             &
     &                      f_ero_nbfrag-f_mass(iv1-1))/(f_mass(iv1)-   &
     &                      f_mass(iv1-1))
                 ENDIF

              ELSEIF ((f_mass(iv2)-f_mass(f_ero_iv)*f_ero_nbfrag) .gt.  &
     &            f_mass(iv1) .and. (f_mass(iv2)-f_mass(f_ero_iv)*      &
     &            f_ero_nbfrag) .lt. f_mass(iv1+1)) THEN

                 f_weight=1.0_sp-(f_mass(iv2)-f_mass(f_ero_iv)*         &
     &                    f_ero_nbfrag-f_mass(iv1))/(f_mass(iv1+1)-     &
     &                    f_mass(iv1))

              ELSE
                 f_weight=0.0_sp
              ENDIF

              ! update for eroded floc masses 

              f_g3(iv2,iv1)=f_g3(iv2,iv1)+f_ero_frac*f_weight*f_beta*   &
     &            f_diam(iv2)*(max(eps,(f_diam(iv2)-f_dp0))/f_dp0)**    &
     &            (3.0_sp-f_nf)*(f_mass(iv2)-f_mass(f_ero_iv)*          &
     &            f_ero_nbfrag)/f_mass(iv1)

              IF (iv1 .eq. f_ero_iv) THEN

                 f_g3(iv2,iv1)=f_g3(iv2,iv1)+f_ero_frac*f_beta*         &
     &              f_diam(iv2)*(max(eps,(f_diam(iv2)-f_dp0))/f_dp0)**  &
     &              (3.0_sp-f_nf)*f_ero_nbfrag*f_mass(f_ero_iv)/        &
     &              f_mass(iv1)
              ENDIF
           ENDIF
        ENDIF ! condition on dfragmax
       ENDDO
      ENDDO

  !********************************************************************************
  !  Shear agregation : LOSS : f_l1
  !********************************************************************************

      DO iv1=1,NCS
       DO iv2=1,NCS

        IF(iv2 .eq. iv1) THEN
           mult=2.0_sp
        ELSE
           mult=1.0_sp
        ENDIF

        f_l1_sh(iv2,iv1)=mult*f_alpha*f_coll_prob_sh(iv2,iv1) 
        f_l1_ds(iv2,iv1)=mult*f_alpha*f_coll_prob_ds(iv2,iv1) 

       ENDDO
      ENDDO

  !********************************************************************************
  !  Shear fragmentation : LOSS : f_l2
  !********************************************************************************

      DO iv1=1,NCS
       f_l3(iv1)=0.0_sp
       IF (f_diam(iv1).gt.dfragmax) THEN
        ! shear fragmentation
           f_l3(iv1)=f_l3(iv1)+(1.0_sp-f_ero_frac)*f_beta*f_diam(iv1)*  &
     &            ((f_diam(iv1)-f_dp0)/f_dp0)**(3.0_sp-f_nf)

        ! shear erosion
        IF ((f_mass(iv1)-f_mass(f_ero_iv)*f_ero_nbfrag) .gt.            &
     &            f_mass(f_ero_iv)) THEN
           f_l3(iv1)=f_l3(iv1)+f_ero_frac*f_beta*f_diam(iv1)*           &
     &            ((f_diam(iv1)-f_dp0)/f_dp0)**(3.0_sp-f_nf)
        ENDIF
       ENDIF
      ENDDO

   IF(MSR)THEN
      WRITE(*,*) ' '
      write(*,*) 'Sum of kernal coefficients:'
      write(*,*) 'f_coll_prob_sh',sum(f_coll_prob_sh)
      write(*,*) 'f_coll_prob_ds',sum(f_coll_prob_ds)
      write(*,*) 'f_g1_sh',sum(f_g1_sh)
      write(*,*) 'f_g1_ds',sum(f_g1_ds)
      write(*,*) 'f_l1_sh',sum(f_l1_sh)
      write(*,*) 'f_l1_ds',sum(f_l1_ds)
      write(*,*) 'f_g3',sum(f_g3)
      write(*,*) 'f_l3',sum(f_l3)
      WRITE(*,*) ' '
      WRITE(*,*) '*** END FLOCMOD INIT *** '   
   END IF

   if(dbg_set(dbg_sbr)) write(ipt,*) "End: initialize_sedflocs_param"

!  END FORMER SUBROUTINE flocmod_kernels

  RETURN

  END SUBROUTINE initialize_sedflocs_param




!=================================================================== ===
!                                                                      !
!  This routine computes floc transformation.                          !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Verney, R., Lafite, R., Claude Brun-Cottan, J., & Le Hir, P. (2011).!
!    Behaviour of a floc population during a tidal cycle: laboratory   !
!    experiments and numerical modelling. Continental Shelf Research,  !
!    31(10), S64-S83.                                                  !
!=======================================================================
!
  Subroutine calc_sed_flocmod

  implicit none

  integer :: i, indx, ised, j, k, ks

!
! Variable declarations for floc model
!
  integer  :: iv1
  real(sp), dimension(0:kb) :: Hz_inv
  real(sp) :: Gval,diss,mneg,dttemp,f_dt
  real(sp) :: dt1,f_csum,epsilon8
  real(sp) :: cvtotmud,tke_av, gls_av, exp1, exp2, exp3, ustr2
  real(sp), dimension(0:kb,ncs) :: susmud
  real(sp), dimension(0:kb,0:mt) :: f_davg
  real(sp), dimension(0:kb,0:mt) :: f_d50
  real(sp), dimension(0:kb,0:mt) :: f_d90
  real(sp), dimension(0:kb,0:mt) :: f_d10
  real(sp), dimension(1:NCS)  :: cv_tmp,NNin,NNout
! f_mneg_param : negative mass tolerated to avoid small sub time step (g/l)
  real(sp), parameter :: f_mneg_param=0.000_sp
  real(sp), parameter :: vonKar = 0.41_sp

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_flocmod "


      epsilon8=epsilon(1.0)
!     epsilon8=1.e-8
!
!--------------------------------------------------------------------------
! * Executable part
!
  do i=1,m 
!
!  Extract mud variables from tracer arrays, place them into
!  scratch arrays, and restrict their values to be positive definite.
!
!  Compute inverse thicknesses of vertical layer to avoid repeated divisions. 
     DO k=1,kbm1
        Hz_inv(k)=1.0_sp/(dt(i)*dz(i,k))
     END DO
     DO ised=1,NCS
        DO k=1,kbm1
           susmud(k,ised)=sed(ised)%conc(i,k)*Hz_inv(k)
         ENDDO
     ENDDO
 
! min concentration below which flocculation processes are not calculated
!      f_clim=0.001_sp 
     exp1 = 3.0_sp+gls_p/gls_n
     exp2 = 1.5_sp+gls_m/gls_n
     exp3 = -1.0_sp/gls_n

     DO k=kbm1,1,-1
 
        f_dt=DTsed
        dttemp=0.0_sp

        ! concentration of all mud classes in one grid cell
        cvtotmud=0.0_sp
        DO ised=1,NCS
           cv_tmp(ised)=susmud(k,ised) 
           cvtotmud=cvtotmud+cv_tmp(ised)

           NNin(ised)=cv_tmp(ised)/sedflocs%f_mass(ised)
        ENDDO

        DO iv1=1,NCS
           IF (NNin(iv1).lt.0.0_sp) THEN
             WRITE(*,*) '***************************************'
             WRITE(*,*) 'CAUTION, negative mass at node i,k :',  &
     &                          i,k
             WRITE(*,*) '***************************************'        
           ENDIF
        ENDDO

        IF (cvtotmud .gt. f_clim) THEN

          if(FLOC_TURB_DISS .and. (.not.FLOC_BBL_DISS))then
!
!ALA dissipation from turbulence clossure
!

!  tke : Turbulent energy squared (m2/s2)                        !
!  gls : Turbulent energy squared times turbulent length scale   !

#  if defined (GOTM)   
            ! use gotm's dissipation rate with caution.
            diss = teps(i,k+1)
#  else
    ! The original turbulence closure scheme is applied and has been tested.
	    IF (k.eq.kbm1) THEN
               tke_av = q2(i,k)
               gls_av = q2l(i,k)

	    ELSEIF (k.eq.1) THEN
               tke_av = q2(i,k)
               gls_av = q2l(i,k)
	    ELSE 
               tke_av = q2(i,k+1)+q2(i,k)
               gls_av = q2l(i,k+1)+q2l(i,k)
            ENDIF
            !tke_av = q2(i,k)
            !gls_av = q2l(i,k)
            diss = gls_cmu0**exp1*tke_av**exp2*gls_av**exp3
#  endif
          elseif(FLOC_BBL_DISS .and. (.not.FLOC_TURB_DISS))then
!
! ALA dissipation from wavecurrent bottom stress
! NOT READY FOR PRIME TIME
! NEEDS VERTICAL DISTRIBUTION
!
              WRITE(*,*) 'CAUTION :'
              WRITE(*,*) 'CHOOSE A DISSIPATION MODEL FOR FLOCS'
              WRITE(*,*) 'BBL MODEL IS NOT AVAILABLE'
              WRITE(*,*) 'SIMULATION STOPPED'
              STOP
          else 
            diss = epsilon8
            IF (l_testcase) THEN
!             if (j.eq.1.and.i.eq.1.and.k.eq.1) then
!               WRITE(*,*) 'VERNEY ET AL TESTCASE FOR FLOCS'
!             endif
            ELSE    
              WRITE(*,*) 'CAUTION :'
              WRITE(*,*) 'CHOOSE A DISSIPATION MODEL FOR FLOCS'
              WRITE(*,*) 'SIMULATION STOPPED'
              STOP
            ENDIF                 
          endif	
      
          CALL flocmod_comp_g(k,i,Gval,diss) 

          DO WHILE (dttemp .le. DTsed)
             CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt)
             CALL flocmod_mass_control(NNout,mneg)
             IF (mneg .gt. f_mneg_param) THEN
                DO WHILE (mneg .gt. f_mneg_param)
                  f_dt=MIN(f_dt/2.0_sp,DTsed-dttemp)
                  IF (f_dt.lt.epsilon8) THEN
	             CALL flocmod_mass_redistribute(NNin)
                     dttemp=DTsed
                     exit
                  ENDIF
                  CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt)
                  CALL flocmod_mass_control(NNout,mneg)   
                ENDDO
             ELSE

                IF (f_dt.lt.DTsed) THEN
                   DO while (mneg .lt.f_mneg_param) 
                      IF (dttemp+f_dt .eq. DTsed) THEN
                         CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt)
                         exit
                      ELSE
                         dt1=f_dt
                         f_dt=MIN(2.0_sp*f_dt,DTsed-dttemp)
                         CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt)
                         CALL flocmod_mass_control(NNout,mneg)
                         IF (mneg .gt. f_mneg_param) THEN 
                           f_dt=dt1
                           CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt)
                           exit
                         ENDIF
                      ENDIF
                   ENDDO
                ENDIF
             ENDIF
             dttemp=dttemp+f_dt

             NNin(:)=NNout(:) ! update new Floc size distribution
             ! redistribute negative masses IF any on positive classes, 
             ! depends on f_mneg_param
	     CALL flocmod_mass_redistribute(NNin)

             IF (dttemp.eq.DTsed) exit

          ENDDO ! loop on full dt

        ENDIF ! only if cvtotmud > f_clim
        ! update mass concentration for all mud classes
        DO ised=1,NCS
           susmud(k,ised)=NNin(ised)*sedflocs%f_mass(ised)
        ENDDO

     ENDDO
!
!-----------------------------------------------------------------------
!  Update global tracer variables.
!-----------------------------------------------------------------------
!
     DO ised=1,NCS
        indx = idsed(ised)
        DO k=1,kbm1
             sed(ised)%conc(i,k)=susmud(k,ised)*(dt(i)*dz(i,k))
        ENDDO
     ENDDO
      
      
  End do
!      WRITE(*,*) 'END flocmod_main'

  if(dbg_set(dbg_sbr)) write(ipt,*) "End:  calc_sed_flocmod"

  End Subroutine calc_sed_flocmod




!!===========================================================================
      SUBROUTINE flocmod_collfrag(Gval,f_g4,f_l4)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_collfrag  ***
  !&E
  !&E ** Purpose : computation of collision fragmentation term, 
  !&E              based on McAnally and Mehta, 2001
  !&E
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_comp_fsd
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------
  !! * Modules used
  !
      implicit none 
 
      real(sp),intent(in) :: Gval

  !! * Local declarations
      integer      :: iv1,iv2,iv3
      real(sp) :: f_fp,f_fy,f_cfcst,gcolfragmin,gcolfragmax
      real(sp) :: gcolfragiv1,gcolfragiv2,f_weight,mult
      real(sp) :: cff1, cff2
      real(sp),DIMENSION(NCS,NCS,NCS) :: f_g4
      real(sp),DIMENSION(NCS,NCS) :: f_l4
      real(sp),DIMENSION(NCS) :: fdiam,fmass
 
      if(dbg_set(dbg_sbr)) write(ipt,*) "Start:  flocmod_collfrag"


! shorten names
      fdiam=SEDFLOCS%f_diam
      fmass=SEDFLOCS%f_mass
!  !!--------------------------------------------------------------------------
  !! * Executable part

      f_fp=0.1_sp
      f_fy=1e-10
      f_cfcst=3.0_sp/16.0_sp
      f_g4(1:NCS,1:NCS,1:NCS)=0.0_sp
      f_l4(1:NCS,1:NCS)=0.0_sp
      cff1=2.0_sp/(3.0_sp-f_nf)
      cff2=1.0_sp/rhoref

      DO iv1=1,NCS
       DO iv2=1,NCS
        DO iv3=iv2,NCS
           ! fragmentation after collision probability based on 
	   ! Gval for particles iv2 and iv3
           ! gcolfrag=(collision induced shear) / (floc strength)

           gcolfragmin=2.0_sp*(Gval*(fdiam(iv2)+fdiam(iv3)))            &
     &          **2.0_sp*fmass(iv2)*fmass(iv3)/(pi*f_fy*f_fp*           &
     &          fdiam(iv3)**2.0_sp*(fmass(iv2)+fmass(iv3))              &
     &          *((SEDFLOCS%f_rho(iv3)-rhoref)*cff2)**cff1)

           gcolfragmax=2.0_sp*(Gval*(fdiam(iv2)+fdiam(iv3)))            &
     &          **2.0_sp*fmass(iv2)*fmass(iv3)/(pi*f_fy*f_fp*           &
     &          fdiam(iv2)**2.0_sp*(fmass(iv2)+fmass(iv3))              &
     &          *((SEDFLOCS%f_rho(iv2)-rhoref)*cff2)**cff1)


           ! first case : iv3 not eroded, iv2 eroded forming 2 particles
	   ! : iv3+f_cfcst*iv2 / iv2-f_cfcst*iv2
           IF (gcolfragmin.lt.1.0_sp .and. gcolfragmax.ge.1.0_sp) THEN  

              IF (((fmass(iv3)+f_cfcst*fmass(iv2)).gt.fmass(iv1-1))     &
     &             .and. ((fmass(iv3)+f_cfcst*fmass(iv2)).le.           &
     &             fmass(iv1))) THEN

                 f_weight=((fmass(iv3)+f_cfcst*fmass(iv2)-              &
     &                     fmass(iv1-1))/(fmass(iv1)-fmass(iv1-1)))

              ELSEIF (fmass(iv3)+f_cfcst*fmass(iv2).gt.fmass(iv1)       &
     &              .and. fmass(iv3)+f_cfcst*fmass(iv2).lt.             &
     &              fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_sp
                 ELSE

                    f_weight=1.0_sp-((fmass(iv3)+f_cfcst*fmass(iv2)-    &
     &                      fmass(iv1))/(fmass(iv1+1)-fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_sp
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &            (SEDFLOCS%f_coll_prob_sh(iv2,iv3))*(fmass(iv3)+   &
     &            f_cfcst*fmass(iv2))/fmass(iv1)

              IF (fmass(iv2)-f_cfcst*fmass(iv2).gt.fmass(iv1-1)         &
     &            .and. fmass(iv2)-f_cfcst*fmass(iv2).le.               &
     &            fmass(iv1)) THEN

                 f_weight=((fmass(iv2)-f_cfcst*fmass(iv2)-              &
     &                     fmass(iv1-1))/(fmass(iv1)-fmass(iv1-1)))

              ELSEIF (fmass(iv2)-f_cfcst*fmass(iv2).gt.fmass(iv1)       &
     &            .and. fmass(iv2)-f_cfcst*fmass(iv2).lt.               &
     &            fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_sp
                 ELSE

                    f_weight=1.0_sp-((fmass(iv2)-f_cfcst*fmass(iv2)-    &
     &                       fmass(iv1))/(fmass(iv1+1)-fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_sp
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &            (SEDFLOCS%f_coll_prob_sh(iv2,iv3))*(fmass(iv2)-   &
     &            f_cfcst*fmass(iv2))/fmass(iv1)


              ! second case : iv3 eroded and iv2 eroded forming 3 particles : 
	      !iv3-f_cfcst*iv3 / iv2-f_cfcst*iv2 / f_cfcst*iv3+f_cfcst*iv2
           ELSEIF (gcolfragmin.ge.1.0_sp .and. gcolfragmax.ge.          &
     &            1.0_sp) THEN  

              IF (f_cfcst*fmass(iv2)+f_cfcst*fmass(iv3).gt.             &
     &            fmass(iv1-1) .and. f_cfcst*fmass(iv2)+f_cfcst*        &
     &            fmass(iv3).le.fmass(iv1)) THEN

                 f_weight=((f_cfcst*fmass(iv2)+f_cfcst*fmass(iv3)-      &
     &                   fmass(iv1-1))/(fmass(iv1)-fmass(iv1-1)))

              ELSEIF (f_cfcst*fmass(iv2)+f_cfcst*fmass(iv3).gt.         &
     &            fmass(iv1) .and. f_cfcst*fmass(iv2)+f_cfcst*          & 
     &            fmass(iv3).lt.fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_sp
                 ELSE

                    f_weight=1.0_sp-((f_cfcst*fmass(iv2)+f_cfcst*       &
     &                 fmass(iv3)-fmass(iv1))/(fmass(iv1+1)-            &
     &                 fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_sp
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &            (SEDFLOCS%f_coll_prob_sh(iv2,iv3))*(f_cfcst*      &
     &            fmass(iv2)+f_cfcst*fmass(iv3))/fmass(iv1)

              IF ((1.0_sp-f_cfcst)*fmass(iv2).gt.fmass(iv1-1)           &
     &            .and. (1.0_sp-f_cfcst)*fmass(iv2).le.                 &
     &            fmass(iv1)) THEN

                 f_weight=((1.0_sp-f_cfcst)*fmass(iv2)-                 &
     &                    fmass(iv1-1))/(fmass(iv1)-fmass(iv1-1))

              ELSEIF ((1.0_sp-f_cfcst)*fmass(iv2).gt.fmass(iv1)         &
     &            .and. (1.0_sp-f_cfcst)*fmass(iv2).lt.                 &
     &            fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_sp
                 ELSE

                    f_weight=1.0_sp-(((1.0_sp-f_cfcst)*fmass(iv2)-      &
     &                      fmass(iv1))/(fmass(iv1+1)-fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_sp
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &                  (SEDFLOCS%f_coll_prob_sh(iv2,iv3))*         &
     &                  ((1.0_sp-f_cfcst)*fmass(iv2))/fmass(iv1) 


              IF ((1.0_sp-f_cfcst)*fmass(iv3).gt.fmass(iv1-1) .and.     &
     &            (1.0_sp-f_cfcst)*fmass(iv3).le.fmass(iv1)) THEN

                 f_weight=((1.0_sp-f_cfcst)*fmass(iv3)-fmass(iv1-1))    &
     &                   /(fmass(iv1)-fmass(iv1-1))

              ELSEIF ((1.0_sp-f_cfcst)*fmass(iv3).gt.fmass(iv1)         &
     &            .and. (1.0_sp-f_cfcst)*fmass(iv3).lt.                 &
     &            fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_sp
                 ELSE

                    f_weight=1.0_sp-(((1.0_sp-f_cfcst)*fmass(iv3)-      &
     &                       fmass(iv1))/(fmass(iv1+1)-fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_sp
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &                  (SEDFLOCS%f_coll_prob_sh(iv2,iv3))*         &
     &                  ((1.0_sp-f_cfcst)*fmass(iv3))/fmass(iv1)


           ENDIF ! end collision test case
        ENDDO
       ENDDO
      ENDDO

      DO iv1=1,NCS
       DO iv2=1,NCS

        gcolfragiv1=2.0_sp*(Gval*(fdiam(iv1)+fdiam(iv2)))**2.0_sp*      &
     &            fmass(iv1)*fmass(iv2)/(pi*f_fy*f_fp*fdiam(iv1)        &
     &            **2.0_sp*(fmass(iv1)+fmass(iv2))*                     &
     &            ((SEDFLOCS%f_rho(iv1)-                            &
     &            rhoref)*cff2)**cff1)

        gcolfragiv2=2.0_sp*(Gval*(fdiam(iv1)+fdiam(iv2)))**2.0_sp*      &
     &            fmass(iv1)*fmass(iv2)/(pi*f_fy*f_fp*fdiam(iv2)        &
     &            **2.0_sp*(fmass(iv1)+fmass(iv2))*                     &
     &            ((SEDFLOCS%f_rho(iv2)-                            &
     &            rhoref)*cff2)**cff1)    

        mult=1.0_sp
        IF (iv1.eq.iv2) mult=2.0_sp

        IF (iv1.eq.MAX(iv1,iv2) .and. gcolfragiv1.ge.1.0_sp) THEN
           f_l4(iv2,iv1)=f_l4(iv2,iv1)+mult*                            &
     &            (SEDFLOCS%f_coll_prob_sh(iv1,iv2))
        ELSEIF (iv1.eq.MIN(iv1,iv2) .and. gcolfragiv2.ge.1.0_sp) THEN
           f_l4(iv2,iv1)=f_l4(iv2,iv1)+mult*                            &
     &            (SEDFLOCS%f_coll_prob_sh(iv1,iv2))
        ENDIF

       ENDDO
      ENDDO

      f_g4(1:NCS,1:NCS,1:NCS)=                                          &
     &            f_g4(1:NCS,1:NCS,1:NCS)*f_collfragparam
      f_l4(1:NCS,1:NCS)=f_l4(1:NCS,1:NCS)*f_collfragparam


      if(dbg_set(dbg_sbr)) write(ipt,*) "End:  flocmod_collfrag"

      RETURN
      END SUBROUTINE flocmod_collfrag

!!===========================================================================
      SUBROUTINE flocmod_comp_fsd(NNin,NNout,Gval,f_dt)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_comp_fsd  ***
  !&E
  !&E ** Purpose : computation of floc size distribution
  !&E
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_main
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------
  !
      implicit none 

  !! * Arguments
      real(sp),intent(in) :: Gval
      real(sp),dimension(1:NCS),intent(in)  :: NNin
      real(sp),dimension(1:NCS),intent(out) :: NNout
      real(sp),intent(in) :: f_dt

  !! * Local declarations
      integer      :: iv1,iv2,iv3
      real(sp) :: tmp_g1,tmp_g3,tmp_l1,tmp_l3,tmp_l4,tmp_g4
      real(sp),dimension(1:NCS,1:NCS,1:NCS)     :: f_g1_tmp,f_g4
      real(sp),dimension(1:NCS,1:NCS)           :: f_l1_tmp,f_l4
      
  !!--------------------------------------------------------------------------
  !! * Executable part


      if(dbg_set(dbg_sbr)) write(ipt,*) "Start: flocmod_comp_fsd "

      tmp_g1=0.0_sp
      tmp_g3=0.0_sp
      tmp_g4=0.0_sp
      tmp_l1=0.0_sp
      tmp_l3=0.0_sp
      tmp_l4=0.0_sp    
      f_g1_tmp(1:NCS,1:NCS,1:NCS)=0.0_sp
      f_l1_tmp(1:NCS,1:NCS)=0.0_sp

      IF (l_COLLFRAG) CALL flocmod_collfrag(Gval,f_g4,f_l4)

      DO iv1=1,NCS
       DO iv2=1,NCS
        DO iv3=1,NCS
           IF (l_ASH) THEN
              f_g1_tmp(iv2,iv3,iv1)=f_g1_tmp(iv2,iv3,iv1)+              &
     &            SEDFLOCS%f_g1_sh(iv2,iv3,iv1)*Gval   
           ENDIF
           IF (l_ADS) THEN
              f_g1_tmp(iv2,iv3,iv1)=f_g1_tmp(iv2,iv3,iv1)+              &
     &            SEDFLOCS%f_g1_ds(iv2,iv3,iv1)   
           ENDIF

           tmp_g1=tmp_g1+(NNin(iv3)*(f_g1_tmp(iv2,iv3,iv1))*NNin(iv2))

           IF (l_COLLFRAG) THEN
              tmp_g4=tmp_g4+(NNin(iv3)*                                 &
     &            (SEDFLOCS%f_g4(iv2,iv3,iv1)*Gval)*NNin(iv2))
           ENDIF
        ENDDO

        tmp_g3=tmp_g3+SEDFLOCS%f_g3(iv2,iv1)*NNin(iv2)*             &
     &            Gval**1.5_sp

        IF (l_ASH) THEN
           f_l1_tmp(iv2,iv1)=f_l1_tmp(iv2,iv1)+                         &
     &            SEDFLOCS%f_l1_sh(iv2,iv1)*Gval
        ENDIF
        IF (l_ADS) THEN
           f_l1_tmp(iv2,iv1)=f_l1_tmp(iv2,iv1)+                         &
     &            SEDFLOCS%f_l1_ds(iv2,iv1)*Gval
        ENDIF

        tmp_l1=tmp_l1+(f_l1_tmp(iv2,iv1))*NNin(iv2)

        IF (l_COLLFRAG) THEN
           tmp_l4=tmp_l4+(SEDFLOCS%f_l4(iv2,iv1)*Gval)*NNin(iv2)
        ENDIF
       ENDDO

       tmp_l1=tmp_l1*NNin(iv1)
       tmp_l4=tmp_l4*NNin(iv1)
       tmp_l3=SEDFLOCS%f_l3(iv1)*Gval**1.5_sp*NNin(iv1)

       NNout(iv1)=NNin(iv1)+ f_dt*(tmp_g1+tmp_g3+tmp_g4-(tmp_l1+        &
     &            tmp_l3+tmp_l4))

       tmp_g1=0.0_sp
       tmp_g3=0.0_sp
       tmp_g4=0.0_sp
       tmp_l1=0.0_sp
       tmp_l3=0.0_sp
       tmp_l4=0.0_sp    
      ENDDO

      if(dbg_set(dbg_sbr)) write(ipt,*) "End: flocmod_comp_fsd "
      RETURN
      END SUBROUTINE flocmod_comp_fsd



!!===========================================================================
      SUBROUTINE flocmod_mass_control(NN,mneg)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_mass_control  ***
  !&E
  !&E ** Purpose : Compute mass in every class after flocculation and 
  !&E              returns negative mass if any
  !&E
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_main
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------

!
      implicit none 
 
  !! * Local declarations
      integer      :: iv1
      real(sp),intent(out)     :: mneg
      real(sp),dimension(1:NCS),intent(in)     :: NN
      !real(sp),DIMENSION(0:NCS+1)      :: f_mass

  !!--------------------------------------------------------------------------
  !! * Executable part
      if(dbg_set(dbg_sbr)) write(ipt,*) "Start: flocmod_mass_control "
      mneg=0.0_sp

      DO iv1=1,NCS
       IF (NN(iv1).lt.0.0_sp) THEN
         mneg=mneg-NN(iv1)*SEDFLOCS%f_mass(iv1)
       ENDIF
      ENDDO

      if(dbg_set(dbg_sbr)) write(ipt,*) "End: flocmod_mass_control "
      RETURN
      END SUBROUTINE flocmod_mass_control

!!===========================================================================
      SUBROUTINE flocmod_mass_redistribute(NN)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_mass_redistribute  ***
  !&E
  !&E ** Purpose : based on a tolerated negative mass parameter, negative masses  
  !&E              are redistributed linearly towards remaining postive masses 
  !&E              and negative masses are set to 0
  !&E                   
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_main
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------
  !
      implicit none 


  !! * Local declarations
      integer      :: iv
      real(sp)     :: npos
      real(sp)     :: mneg
      real(sp),dimension(1:NCS),intent(inout)     :: NN
      real(sp),dimension(1:NCS)                   :: NNtmp
      !real(sp),DIMENSION(0:NCS+1)      :: f_mass
  !!--------------------------------------------------------------------------
  !! * Executable part

      if(dbg_set(dbg_sbr)) write(ipt,*) "Start:  flocmod_mass_redistribute"

      mneg=0.0_sp
      npos=0.0_sp
      NNtmp(:)=NN(:)

      DO iv=1,NCS
       IF (NN(iv).lt.0.0_sp) THEN
         mneg=mneg-NN(iv)*SEDFLOCS%f_mass(iv)
         NNtmp(iv)=0.0_sp
       ELSE
         npos=npos+1.0_sp
       ENDIF
      ENDDO

      IF (mneg.gt.0.0_sp) THEN 
       IF (npos.eq.0.0_sp) THEN
         WRITE(*,*) 'CAUTION : all floc sizes have negative mass!'
         WRITE(*,*) 'SIMULATION STOPPED'
         STOP    
       ELSE
         DO iv=1,NCS
           IF (NN(iv).gt.0.0_sp) THEN
              NN(iv)=NN(iv)-mneg/sum(NNtmp)*NN(iv)/                     &
     &            SEDFLOCS%f_mass(iv)
           ELSE
              NN(iv)=0.0_sp
           ENDIF

         ENDDO

       ENDIF
      ENDIF

      if(dbg_set(dbg_sbr)) write(ipt,*) "End:  flocmod_mass_redistribute"
      RETURN
      END SUBROUTINE flocmod_mass_redistribute

!!===========================================================================    
  SUBROUTINE flocmod_comp_g(k,i,Gval,diss)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_comp_g  ***
  !&E
  !&E ** Purpose : compute shear rate to estimate shear aggregation and erosion  
  !&E 
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_main
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------
  !
  implicit none 

  !! * Local declarations
  integer,  intent(in)      :: k,i
  real(sp),intent(out)     :: Gval
  real(sp)     :: htn,ustar,z,diss,nueau
! l_testcase - if .TRUE. sets G(t) to values from lab experiment
!      logical, parameter :: l_testcase = .TRUE.
  !!--------------------------------------------------------------------------
  !! * Executable part

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: flocmod_comp_g "

  !    nueau=1.06e-6_sp
  !    ustar=sqrt(tenfon(i)/rhoref)
  !    htn=h0(i)+ssh(i)
  !    z=(1.0_sp+sig(k))*htn

!
! ALA from CRS
!
  nueau = 1.5E-6_sp
  IF (l_testcase) THEN
    ! reproducing flocculation experiment Verney et al., 2011
!    Gval=0.0_sp
!    IF (time .lt. 7201.0_sp) THEN
!      Gval=1.0_sp
!    ELSEIF (time .lt. 8401.0_sp) THEN
!      Gval=2.0_sp
!    ELSEIF (time .lt. 9601.0_sp) THEN
!      Gval=3.0_sp  
!    ELSEIF (time .lt. 10801.0_sp) THEN
!      Gval=4.0_sp
!    ELSEIF (time .lt. 12601.0_sp) THEN
!      Gval=12.0_sp
!    ELSEIF (time .lt. 13801.0_sp) THEN
!      Gval=4.0_sp
!    ELSEIF (time .lt. 15001.0_sp) THEN
!      Gval=3.0_sp
!    ELSEIF (time .lt. 16201.0_sp) THEN
!      Gval=2.0_sp
!    ELSEIF (time .lt. 21601.0_sp) THEN
!      Gval=1.0_sp
!    ELSEIF (time .lt. 25201.0_sp) THEN
!      Gval=0.0_sp
!    ELSEIF (time .lt. 30601.0_sp) THEN
!      Gval=1.0_sp
!    ELSEIF (time .lt. 31801.0_sp) THEN
!      Gval=2.0_sp                     
!    ELSEIF (time .lt. 33001.0_sp) THEN
!      Gval=3.0_sp       
!    ELSEIF (time .lt. 34201.0_sp) THEN
!      Gval=4.0_sp
!    ELSEIF (time .lt. 36001.0_sp) THEN
!      Gval=12.0_sp
!    ELSEIF (time .lt. 37201.0_sp) THEN
!      Gval=4.0_sp
!    ELSEIF (time .lt. 38401.0_sp) THEN
!      Gval=3.0_sp
!    ELSEIF (time .lt. 39601.0_sp) THEN
!      Gval=2.0_sp
!    ELSEIF (time .lt. 45001.0_sp) THEN
!      Gval=1.0_sp
!    ELSEIF (time .lt. 48601.0_sp) THEN
!      Gval=0.0_sp
!    ELSEIF (time .lt. 54001.0_sp) THEN
!      Gval=1.0_sp
!    ELSEIF (time .lt. 55201.0_sp) THEN
!      Gval=2.0_sp                     
!    ELSEIF (time .lt. 56401.0_sp) THEN
!      Gval=3.0_sp       
!    ELSEIF (time .lt. 57601.0_sp) THEN
!      Gval=4.0_sp
!    ELSE 
!      Gval=12.0_sp
!    ENDIF
  ELSE
    Gval=sqrt(diss/nueau)
  ENDIF
! NO KLUDGE
!      Gval = 12.0_sp

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: flocmod_comp_g "
  RETURN

  END SUBROUTINE flocmod_comp_g
      
      
# endif

End Module Mod_FlocMod

